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Abstract 

We have implemented the log-conformation method for two-dimensional viscoelastic flow in COMSOL, a 
commercial high-level finite element package. The code is verified for an Oldroyd-B fluid flowing past a 
confined cylinder. We are also able to describe the well-known bistability of the viscoelastic flow in a cross¬ 
slot geometry for a FENE-CR fluid, and we describe the changes required for performing simulations with 
the Phan-Thien-Tanner (PTT), Giesekus and FENE-P models. Finally, we calculate the flow of a FENE- 
CR fluid in a geometry with three in- and outlets. The implementation is included in the supplementary 
material, and we hope that it can inspire new as well as experienced researchers in the field of differential 
constitutive equations for viscoelastic flow. 

1 Introduction 

The flow of viscoelastic fluids in complex geometries has many practical applications due to the fact that 
any dissolved elastic component, be it polymers or biological molecules, will give rise to viscoelastic effects. 
Such effects require the use of models that take the fluid’s memory of past deformations into account. 
Differential constitutive equations is able to do just that, and they provide a good quantitative agreement 
with experiments pQ. Complex geometries call for the use of numerical methods, and therefore significant 
effort has been devoted towards various discretization techniques as well as model reformulations. The 
appearance of singularities, that caused codes to break down in smooth geometries at moderate elasticity, 
led to the definition of a ’’High Weissenberg Number Problem” (HWNP). This was effectively solved, when 
loss of positive definiteness for the conformation tensor was recognized as the cause of the singularity and 
the log-conformation method was introduced as a remedy [2]. 
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Commercial numerical tools provide a simple workflow due to the integration of geometry description, 
unstructured mesh generation, discretization, solvers, and post-processing, which makes the treatment of 
complex geometries effortless compared to research grade code, but to our knowledge there does not exist 
any commercial tool that implements the log-conformation method. This is perhaps due to the fact that 
the method involves the calculation of eigenvectors and eigenvalues for the conformation tensor, which 
complicates the formulation of the constitutive equation. In two (and three) dimensions explicit expressions 
exists for these quantities, and it is thus possible to formulate the governing equations solely in terms of 
functions recognized by commercial simulation packages. We have included a COMSOL implementation of 
the Oldroyd-B model with log-conformation reformulation in the supplementary material, and we hope this 
can help other researchers working with differential constitutive equations. 

The article is structured in four parts: First we discuss implementation details of the log-conformation 
method for the most simple constitutive equations. We then calculate the flow of an Oldroyd-B fluid past a 
confined cylinder using COMSOL for comparison with other works [2]. As a third point we demonstrate the 
ability to predict bistable behavior of a FENE-CR fluid in a cross-slot geometry [3j. Finally we consider a 
geometry with three in- and outlets, which has not previously been described in the context of viscoelastic 
flow, at least to our knowledge. 

2 Log-Conformation Implementation 

In the so-called dumbbell models, the elastic part of a viscoelastic fluid is approximated by two spring- 
connected point masses - an elastic dumbbell. The orientation and elongation of this dumbbell is described 
using the end-to-end vector a with the equilibrium length a eq , see figure fTl The conformation tensor, A, is 



Figure 1: An elastic dumbbell is illustrated with its two spring-connected point masses, and the end-to-end 
vector describing orientation and extension. 


the statistical average, (...) of the dyadic product between the end-to-end vector and itself normalized with 
the squared equilibrium length, 


A = 


(a ( 8 ) a) 


(i) 


The equation system for the Oldroyd-B model is written below. We neglect the effect of temperature 
and compressibility, but keep inertia for the sake of generality. 
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where D is the material derivative, p is the density, v is the velocity vector, t is time, p is the pressure p s is 
the solvent viscosity, r] p is the elastic viscosity, A is the relaxation time, r is the solvent stress tensor and 
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is the elastic stress tensor. The Navier-Stokes equation Q guarantees that Newton’s 2nd law is obeyed, 
while the continuity equation © insures mass conservation. The expression for the elastic stress tensor 
in equation Q is best understood in light of the conformation tensor definition ([!]) with rj p /X acting as a 
spring constant. The evolution of the conformation tensor is described in equation © with relaxation on 
the first line balanced by the upper convected time derivative, which consists of convection in terms of the 
total derivative together with rotation and extension of the conformation tensor as imposed by the velocity 
gradient in square brackets. 

When the velocity vector is approximated by C° continuous polynomials, the velocity gradient, V v 
becomes discontinous. This can cause numerical difficulties [4], and therefore a C° continuous approximation 

G is usually constructed for use in equation (5). Furthermore the zero, r] p v + (V v) T — G — G T ^ is 

added on the right-hand side of equation © to preserve the elliptic nature of the equation for small solvent 
viscosities in what is called the discrete elastic viscous stress splitting (DEVSS). Similarly, we use streamline 
upwind diffusion (SUPG) to introduce an elliptic component in the constitutive equation ([ 5 ]), see section 2.2 


The conformation tensor is symmetric, and the point of the log-conformation method is to exploit this 
property to rewrite equation © to the form 

Ds 

n ~ n( 8 '®=°' 

where s is the matrix logarithm of the conformation tensor. They are related by the following expression 


A = E e= a T 

_ ^lnx ^ln y 

^ln y V\ nx 

g [ e Al 0 

e “ “ [ 0 e Aa 

^11 = I’LT 1 + 

A12 = VlnxVlny (e Al - e A2 ) 

A 22 = vj nx e X2 + v 2 lny e Xl 

where Ai and A 2 are the eigenvalues of s, while v\ nx and v\ ny are the normalized eigenvector components 
belonging to Ai. The orthonormal rotation matrix, R holds the components of the eigenvectors, but we 
exploit that the second eigenvector is just the transpose of the first. Denoting the components of s by sn, 
812 and 822 , the eigenvalues becom^] 


Ai 

A 2 


\ ^22 + sn - \j{s 22 - 8 n ) 2 + 48 f 2 ^ and 
\ ^ 5 22 + S 11 + \J {S22 ~ Sll) 2 + 48f 2 ^j , 


and the components of the first eigenvector are 


vi 


A| — 822 
<§12 


The eigenvector is normalized 


V lx / ^ v l*+ v lj, 

_ «1 V / ^/ v lx+ v l» _ 

The special case 812 = 0 is handled by using R = I, when 1 812 1 < e with e = 10 -12 

1 There also exist explicit expressions for the three-dimensional case. 
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In the following we will express the reaction term II as R • Q • R T as described in [5]. The continuous 
velocity gradient is rewritten in the principal frame 

£ = E T £ E 

Gn = vJ na ,Gii + Vi nx Viyx{Gi2 + G 21 ) + ^i ny G22 

G 12 = v lnccG21 + Vin X Vly X {G22 ~ Gn) — ^i ny Gl2 

G 21 = V^ n;E Gi2 + Vi ncc Vi y;E (G22 — Gn) — V^ ny G21 

G22 = v lna;G22 — Vi 7lIC Vi 2/IC (Gi2 — G21) + V^ n2/ Gn. 

Note that this matrix is not symmetric, and that we have chosen G 12 to refer to the derivative of the velocity 
vector x-component with respect to the //-coordinate. The diagonal components fin and ^22 are 

1 — e~ Xl 

fin = 2Gn--- and 

A 

1 — e~ X2 

^22 = 2G 22 - j -• 

These expressions are specific to the Oldroyd-B model, but it is straight forward to write the expressions for 
the Giesekus, Phan-Thien-Tanner, FENE-P or FENE-CR models instead [5]. 


Giesekus : 

fill — 2Gn — /g ie.(Ai) 
fi 2 2 = 2G22 — /g ie.(A2) 

(6) 

where 

/Gie.(*)= 1 + aGie A (eI_1) ( e * 1) 


PTT(exp.) : 

fin = 2Gu fe PT T (A 1 ,A 2 ) 

n 22 = 2G 22 fepTT ( A Ai ’ A2 ) ’ 


where 

fcpT T (Ai,A 2 ) = e e PTT( A i+A 2 - 3 ) 

(7) 

FENE - P : 
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where 

fc(Ai,A 2 ) ^ . \ \ / 2 

1 — (Ai + A 2 )/max^ 


FENE - CR : 

fin = 2G 11 «fe(A 1 ,A 2 )^ 

fi 22 = 2G 22 fc(A 1; A 2 ) eA A _1 ’ 


where 

^(Ai, A2) 1 /■» | \ \ / 2 

1 — (Ai + A 2 )/< ax 



In the case of the FENE-P and FENE-CR models, the a max is a maximum extension, i.e. Tr(A) < a^ ax . 
The off-diagonal component, fii 2 , equals 

q i2 _ ( ( e M-e^2) ( eAl Gi2 + e A2 G 2 i^ , |Ai - A 2 | > e 
\ G12 + G21 , |Ai — A2I < e 

where we have reused the numerical parameter for the eigenvectors to handle the special case of identical 
eigenvalues. We are now ready to compute the reaction term 

hill = 'Uina^ll — 2Vi nx Vi ny Q,i2 + ^lny^22 
hfl2 = (v\ nx ~ Vlny) fil2 + Vi nx Vi ny (Qn — Q 22 ) 

1122 = ^lncc^22 + 2^i na ,^i ny fii2 + vl ny Qn 

To ease modification of the implementation, we have opted to construct functions for each of the above 
expressions instead of writing out II in terms of expressions involving s and G only. 
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With respect to visualization it is worth noting that s and A are positive definite tensors, and one 
can thus think of them as ’’fields of ellipses”. In other words, the conformation tensor has an extension in 
every point given by its eigenvectors and eigenvalues. The trace of the conformation tensor is the sum of 
the extensions in the principal directions given by the eigenvectors, while the determinant is the product 
of these extensions. The trace and determinant equal the sum and product of the eigenvalues, respectively, 
and they are invariant with respect to rotation of the coordinate system, which make them interesting in 
the context of post-processing/visualizatior0 The last invariant is the angle(s) related to the eigenvectors, 
but the eigenvector associated with the largest eigenvalue is normally oriented with the streamlines, so this 
angle provides little insight. Furthermore the trace and determinant often look very similar, at least on a 
logarithmic scale, and therefore we just plot the trace. 


2.1 Non-Dimensional Equations 

We define the dimensionless spatial coordinates, velocity, time, pressure and stress (x, v, t , p and r g ), 


X = T char x, V 

~ , Lchgj- ~ 

— ^charV, t — £, 

^char 

V 

^char(^7s H" Vp) ~ 

T P 

-^char 

and r 

—e 

/ \ 't’char ~ 

= ( Vs+Vp) T T e , 

-^char 


which give rise to the following dimensionless constants 

P^char-^char ~x~x t \ ^char 

-, We = A-- and 

V T c h ar 

Vs 

Vp T Vs 

The Weissenberg number, We, describes the magnitude of elastic to viscous effects. When either the Weis- 
senberg number goes to zero or the solvent to total viscosity ratio, /3 goes to unity, the model approaches that 
of a Newtonian fluid. Contrarily, the value /3 = 0 corresponds to the upper convected Maxwell model, which 
is characterized by particularly strong viscoelastic effects. Note that we prefer to define a characteristic 
pressure in terms of a characteristic velocity rather than the other way around due the fact that we intend 
to impose fixed flow rates rather than fixed driving pressures (as in 0). 

We have implemented the above approach to the log-conformation formulation for two dimensions in 
COMSOL Multiphysics version 4.3a, and the online version of the manuscript includes a MATLAB script 
for computation of the flow of an Oldroyd-B fluid past a confined cylinder using the COMSOL LiveLink 
interface. We use default settings for everything, except we force update of the Jacobian every time a time 
step is saved. 


Re = 

P = 


2.2 The Oldroyd-B model in Weak Form 


Converting the strong form of the Oldroyd-B model (2l5) to weak form is trivial for the equations that 
relate to the elastic stress Q and continuity equation (|3| as one just multiplies with test functions for the 
elastic stress and pressure, respectively. The evolution equation for the conformation tensor © requires 

(2h me sh/^char)V ' VA^^ 


SUPG stabilization, which means that the equation is multipliecj^j with A^ ( 


rather than A^^, where h mes h is a characteristic size for the mesh. In other words a small amount of 
stream-wise diffusion is added to allow computation, but since the amount of diffusion is scaled with the 
mesh, the implementation still converges towards the solution without any artificial diffusion. Finally there 
is the momentum part of the Navier-Stokes equation which is multiplied with a test function for the 


2 Note that det(A) =a Fh eXi — e^ lXi = e Tr ^, where A i are the eigenvalues of s. 

3 Note that the test functions for elastic stress and conformation tensor are multiplied element wise (Hadamard product), 
denoted with o. 
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velocity and integrated by parts: 


0 = 
0 = 


+ 



’ v test^^ 



° V Vtest dfl 


Dv 

Vf — Re—~ — Da 
= e Dt 



/ r • n • vtest dl 

Jdn 


• Vtest dQ 


( 8 ) 


(9) 


The point is that the solvent stress is integrated by parts (Divergence theorem), but the elastic stress is not. 
This approach has been followed by other researchers so we have adopted it as well. The advantage is that 
the specification of r on open boundaries is slightly more simple than r + r , but in fact we have tested 
both without finding any difference. At outlets we thus usually impose normal velocity, zero pressure and 


T = —Ip. 

=s — 

This can also be used for pressure driven inlets, but in that case one should supplement the constraint on 
the velocity with one on the conformation tensor. We impose the no-slip(v = 0) at walls, which causes 
COMSOL to automatically enforce v tes t = 0? such that the boundary integral in equation vanishes. 


3 Confined Cylinder 

The confined cylinder geometry is shown in figure [2] and it is a popular benchmark geometry due to the 
absence of geometric singularities. It can be used to verify an implementation by comparing the drag over 
the cylinder with a reference. We exploit the symmetry of the problem by modeling only half of the domain. 



Figure 2: The confined cylinder geometry with a (dashed) symmetry axis. The fluid comes in from the left, 
flows past the obstacle and exits through the boundary on the right (t is the unit tangent vector). 


The cylinder radius R, and average inlet velocity are used as characteristic length and velocity, respec¬ 
tively. The blocking ratio is set at 2 (H = 2R), and the outlet length, Li, is fixed at 10 R. In agreement 
with references [2j the solvent to total viscosity ratio is taken as 0.59, and inertia is neglected (Re = 0). The 
dimensionless drag is computed as 


Drag = 

—2 / r h • x T d/, 

J cylinder 

with 

r = 

—pi + Vs Vv + (Vv) T 

+ r 


where n is the outward pointing normal vector. 

We impose the no-slip boundary condition v = 0 at the walls, and fully developed boundary conditions 


6 









at the inlet 


Vinlet 

= f(l-(y/ 2 ) 2 )x 



-11,inlet 

- 1+ ’(VS = 

1 + |f/ 2 We 2 


-12,inlet 

= We^; x) = —fyWe 
dy 4 

and 

( 10 ) 

-22,inlet 

= 1 




Note that the intermediate expressions can be used in the case of a pressure driven setup. At the outlet we 
impose zero pressure together with zero normal viscous stress. 


Poutlet 


r , • n 

=s, out let 


0 and 
(—Ip) • n = 0 


At the symmetry line we impose zero normal velocity. 


v sym ■ n = 0 


( 11 ) 


( 12 ) 


3.1 Solution Details 


In order to find a steady solution one can go through the following steps. 

#1 Solve for the case of a Newtonian fluid. 

#2 Use #1 as initial condition for a transient simulation running for 10 relaxation times. 
#3 Use #2 as initial guess for a non-linear solver to find a steady solution. 


We however choose to vary the Weissenberg number slowly using a regularized step function, st, such that 
a quasi steady state is achieved. 


We = 

lVe s t a rt + (^Ve enc j We 

start) St(t), 

where 


f ° 


t ^ ^start 

st(t) = 

{ 0.5 + 1.54-2 (tf 

5 ^start ^ 

t ^ ^end 


l 1 

5 ^end ^ 

t 

t = 

t (tstart + 4 e nd)/2 



^end ^start 




This way we compute a quasi-steady solution to a range of Weissenberg numbers using a single transient 
simulation. COMSOL uses a fully implicit transient scheme with adaptive stepsize, so slow variation of the 
Weissenberg number does not increase the total computation time (?]. 


3.2 Results 

Since the time over which the Weissenberg number is varied is a numerical parameter, we choose to scale it 
with the mesh size, so convergence towards the steady case is achieved, i.e. 

Tstep = tend ~ tstart = 8000(0.07//l m esh) 

Figure [3] shows the drag as a function of the Weissenberg number for three different discretizations indicating 
convergence up to We = 0.6. We observe good agreement with the reference drag although the simplicity 
of the chosen mesh gives rise to a high computational cost. It is well known, that convergence at We = 0.7 
is difficult, even with anisotropic mesh adaption |8 . In fact experimental [9] as well as theoretical evidence 
m exist to suggest, that a transition to a three dimensional flow pattern occurs, in other words the two 
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Figure 3: The drag on a confined cylinder is plotted as a function of the Weissenberg number for three 
discretizations together with the data of another study [2]. 


We=0.51 



We=0.70 



iioo 
10 Tr(A) 

—h 

Bioo 
10 Tr(A) 


Figure 4: The trace of the conformation tensor is plotted on a logarithmic scale for the confined cylinder in 
the case of We = 0.51 and We = 0.7 (183k DOF). The acceleration of the flow in the wake of the cylinder 
generates a strand of highly elongated dumbbells. 


dimensional symmetric flow becomes unstable. The three-dimensional solution has the flow alternating 
between going above and below the cylinder, and one can think of this solution as having a weaker extensional 
character in the wake of the cylinder. 

When a viscoelastic fluid enters a region, where the extension rate times the relaxation time reaches 
a critical value, the dumbbell extension starts to grow exponentially in time. This kind of non-linearity 
explains the significant difference between the wake of extension at We = 0.51 and We = 0.7 illustrated in 
figure [4] 


4 Bistable Cross-Slot 

The FENE-CR model differs from the Oldroyd-B model in the sense, that A is replaced by A [l — Trace(A) /aj 2 liax ]. 
In other words the relaxation time decreases as the dumbbell extension increases. This puts an upper bound 
a max? on th e t race of the conformation tensor corresponding to a maximum extension. The modification 
prevents an unbounded extension in free stagnation points such as the center of the cross-slot geometry 
shown in figure [5j This geometry is known to exhibit bistability experimentally m as well as numerically 

m- 

In this case we keep the average inlet velocity as characteristic velocity and the channel width H as 
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Figure 5: The cross-slot geometry is shown with the fluid coming in from the left and right, diverging at 
the center, before exiting through to the upper and lower outlets leaving a free stagnation point (red) in the 
center of the geometry. 


characteristic length scale. We impose the no-slip boundary condition at walls, but we have been unable 
to initialize the transient simulation with the correct expressions for the conformation tensor at the inlets. 
Therefore we have used s = el at the inlets rather than values based on the expressions listed in Appendix 
I. Furthermore we had to change the damping strategy of the nonlinear method from the default (constant) 
to automatic. 

At the outlet we impose the same boundary conditions as for the cylinder. We compute the dissi¬ 
pation as an integral over the entire computation domain, 


0 = 



Vv + (Vv) T 


d£l. 


(14) 


where : is the Frobenius product. We choose Re = 0, a^ ax = 100 and /3 = 0.2 and the same transient 
solution procedure as described in section |3.1| We define the center of the domain, 


fic = X G |x|oo < 0.5, 


as we find that the vorticity integrated over this is a good asymmetry parameter. We thus plot the square 
of this versus the Weissenberg number in figure [6^a). The system becomes bistable at We = O.fQ which 
gives rise to a kink in the dissipation as shown in figure [6jb). The solution before the kink [figure UN 
and (d)] has vertical, horizontal as well as 180-degree rotational symmetry, while the solutions after the 
kink 0b), (c), (e) and (f)] only have 180-degree rotational symmetry, and we will refer to these solutions 
with reduced symmetry as ”asymmetric solutions”. The phenomenon can be understood by considering the 
velocity magnitude and conformation tensor trace as plotted in figures [Tja-c) and (d-e) respectively. The 
dumbbells are extended along the vertical axis, which causes a damping that delays the merging of the 
flows, and this gives rise to extra shear and thus also dissipation. When the solution becomes asymmetric, 
the damping is moved to the side of the channel, which lowers the dissipation. One can also think of the 
asymmetry as a flow pattern with less extension in the stagnation point, similar to the three-dimensional 
flow for the confined cylinder. 

4 Based on personal correspondence we have concluded that the results of [5] are based on the FENE-MCR model, so no 
reference exists for the FENE-CR model. 
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(a) 



Figure 6: The square of the vorticity integrated over the center of the domain Q c is plotted as function 
of the Weissenberg number in (a). In (b) it is the dissipation that is plotted, which shows that the point 
of bistability occurs (We = 0.5) at the point of maximum dissipation and thus also maximum hydraulic 
resistance. 

5 Three Inlets 

In this section we introduce the geometry shown in figure [8j It consists of a hexagon with in- and outlets 
attached in an alternating fashion. The geometry always has a vertical symmetry axis, but depending on 
whether a certain angle is equal to tt/3 or not, it can also have 120-degree rotational symmetry. We define 
the distance to the center as in- and outlet lengths, and set equal inlet flow rates as well as outlet pressures. 
We set Lin = 6H and L out = 8L, /3 = 0.1 and s = el at the inlets, but keep all other parameters identical 
to those used for the cross-slot. That is, the channel width and average inlet velocity are still used as 
characteristic length and velocity, respectively. 

If the angle is equal to n/3 radians, a single stagnation point exists in the center of the geometry regardless 
of whether we consider Stokes flow or just the regime of low elasticity. Two stagnation points however appear, 
if the angle is different from n/3 radians. Figure [9] shows the solution for angles of 7r/3, 7 t/ 3.5 and 7r/4 radians 
at We = 0.37 (low) and We = 0.66 (moderate). For an angle of tt/3 radians, the regime of low elasticity (a) 
gives a solution with 120-degree rotational- and vertical symmetry, while for moderate elasticity a solution 
with only 120-degree rotation symmetry appears (b). This latter solution involves the clockwise flow, but by 
symmetry we can argue that there must be an identical solution with counter-clockwise flow. The symmetry 
reduction also occurs for an angle of 7r/3.5 (c-d), in the sense that we go from only vertical symmetry to 
no symmetry. Regardless of the Weissenberg number, there are two stagnations points even though the 
flow rate between them is small at moderate elasticity (d). In other words turning up the elasticity causes 
the flow rate between the stagnation points to decrease significantly. As shown in (e-f) this is different for 
an angle of 7r/4 radians, as the elasticity tends to move the stagnation points apart and increase the flow 
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Figure 7: The magnitude of the velocity vector is plotted for three Weissenberg numbers in the cross geometry 
in (a-c), while red stream lines are plotted on top of the trace of the conformation tensor on a logarithmic 
scale in (d-f). Note that we plot streamlines as contours of the stream function instead of integrating the 
velocity field. 


p=o 



Figure 8: The geometry with three in- and outlets is illustrated with the different length scales and angles 
involved. For an angle of 7r/3 the geometry has a symmetry axis through every channel as well as 120-degree 
rotational symmetry, while an angle different from tt/3 breaks all symmetries except the one around the 
vertical axis. 


rate between them in what we call the ”straight flow solution”, that is the solution has identical symmetry 
properties at low and moderate elasticity. 

Using artificial initial conditions as listed in Appendix II, one can show that the clockwise, counter¬ 
clockwise and the straight flow all exist as stable solutions at moderate elasticity for angles between 7r/4 and 
7r/3.5. This means, that it is the connectivity with the solution at low elasticity that changes. Furthermore 
there is an unstable solution with a vertical symmetry axis for the angles 7r/3 and 7r/3.5. We calculate the 
integral of the vorticity around the center 

U 0 = x G ||x || 2 < 0.5, 
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angle=7i/3, We=0.37 angle=7i/3, We=0.66 



Tr(A) 


( c ) angle=;i/3.5, We=0.37 (d) angle=7t/3.5, We=0.66 



Figure 9: The trace of the conformation tensor is plotted on a logarithmic scale together with streamlines 
(red) for two Weissenberg numbers in a geometry with three inlets. The case of an angle equal to 7t/3 (a-b), 
7t/3.5 (c-d) and 7t/4 (d-e) radians is illustrated (see figure [8]). Although not apparent from the streamlines, 
the upper outlet receives fluid from all three inlets in (d), and the flow between the stagnation points is 
reduced with a factor of 4.5 in (d) compared to (c). Any asymmetry around the vertical axis in (a), (c) and 
(e) is due to the mesh. 


as a measure of asymmetry and plot this as a function of the Weissenberg number for five angles in figure 
[To] (a), while the dissipation is plotted in figure 10 (b). It is clear that the point of maximum dissipation 
coincides with the point of bistability for the three largest angles (7 t/3 , 7 t/3 .5, 7 t/3 .75), as it was also the 
case for the cross-slot in section [4] The dissipation for the two smaller angles also have maxima, but we do 
not believe they play the same role as the maxima for the smaller angles. This is due to the fact that the 
Weissenberg number at which the maximum occurs seems to decreases for both the rotation and straight 
flow solutions. This means that there does not exists an angle for which the two kind of maxima in figure 
[To] (b) meet, although there must exist an angle where all three solutions are connected to the solution at 
low elasticity. 


6 Conclusion 

We have described the viscoelastic Oldroyd-B model with log-conformation reformulation in detail, and 
included a COMSOL implementation in the supplementary material. We have showed, that the code agrees 
with a reference in a benchmark geometry. Furthermore we have demonstrated that the code can be modified 
to describe the flow of a FENE-CR fluid, and that it can reproduce bistability in the cross-slot geometry. 
Finally we computed the flow in a geometry with three in- and outlets, which appears to feature three stable 
solutions. The connectivity of the solutions with the solution at low elasticity seems to depend on an angle 
in the geometry. 
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(a) 



Figure 10: The square of the vorticity integrated over the center of the domain Q 0 is plotted as functions 
of the Weissenberg number for five different angles in (a), while the dissipation is shown (b). The point of 
bistability coincides with the point of maximum dissipation for the three larger angles, as it was also the case 
for the cross-slot in figure [6j The maximum dissipation of the two smaller angles is indicated with crosses. 


References 

[1] F. Baaijens, S.H.A. Selen, H.P.W. Baaijens, G.W.M. Peters, and H.E.H. Meijer. Viscoelastic flow past 
a confined cylinder of a low density polyethylene melt. Journal of non-newtonian fluid mechanics , 68 
(2-3):173-203, 1997. 

[2] M. A. Hulsen, R. Fattal, and R. Kupferman. Flow of viscoelastic fluids past a cylinder at high weissenberg 
number: stabilized simulations using matrix logarithms. Journal of Non-Newtonian Fluid Mechanics , 
127(1):27-39, 2005. 

[3] G.N. Rocha, R.J. Poole, M.A. Alves, and P.J. Oliveira. On extensibility effects in the cross-slot flow 
bifurcation. Journal of Non-Newtonian Fluid Mechanics , 156(l):58-69, 2009. 

[4] A. Fortin, R. Guenette, and R. Pierre. On the discrete evss method. Computer Methods in Applied 
Mechanics and Engineering , 189(1): 121-139, 2000. 

[5] A. Kane, R. Guenette, and A. Fortin. A comparison of four implementations of the log-conformation 
formulation for viscoelastic fluid flows. Journal of Non-Newtonian Fluid Mechanics , 164(1-3):45-50, 
2009. ISSN 0377-0257. 

[6] K. Ejlebjerg Jensen, P. Szabo, and F. Okkels. Topology optimization of viscoelastic rectifiers. Applied 
Physics Letters , 100(23):234102-234102, 2012. 


13 



















[7] COMSOL Multiphysics. 4.3a user’s guide, 2013. 

[8] R. Guenette, A. Fortin, A. Kane, and J.F. Hetu. An adaptive remeshing strategy for viscoelastic fluid 
flow simulations. Journal of Non-Newtonian Fluid Mechanics , 153( 1):34—45, 2008. 

[9] G.H. McKinley, R.C. Armstrong, and R.A. Brown. The wake instability in viscoelastic flow past confined 
circular cylinders. Philosophical Transactions: Physical Sciences and Engineering , 344(1671):265-304, 
1993. 

[10] M. Sahin. Parallel large-scale simulation of viscoelastic fluid flow instabilities. In 17th International 
Workshop on Numerical Methods for Non-Newtonian Flows , 2012. 

[11] P.E. Arratia, CC Thomas, J. Diorio, and JP Gollub. Elastic instabilities of polymer solutions in cross¬ 
channel flow. Physical review letters , 96(14): 144502, 2006. 

[12] RJ Poole, MA Alves, and PJ Oliveira. Purely elastic flow asymmetries. Physical review letters , 99(16): 
164503, 2007. 


Appendix I: FENE-CR Inlet Conditions 

The developed flow conditions for a FENE-CR fluid [^] is listed in the following equations 

Vinlet = -f (1 - (£/2) 2 ) n 


^-11, inlet 

A 12 ,inlet 
^22,inlet 

X 


z -t 

^max 


(«L - y a max + 8X 2 «ax- 2 )) 


= x 1- 


4x 2 

^-11,inlet T 1 


= We 


1, where 
d (v • x) 


dy 


- 3We(i//4)(n • x) 


Appendix II: Artificial Initial Conditions for the Geometry with 
Three In- and Outlets 


One can find stable solutions not connected to the solution at low elasticity using artificial initial conditions. 
A rotating solution for an angle of 7 t/4 (see figure [8]) can in example be found by imposing a rotating force 
in the center of the geometry, 


f {xy-yx)/2 , ||r|| < 0.5/ sin(7r/6) 

\ 0 , jjfjj > 0.5/sin(7r/6) ’ 


where f is the position vector and 11 • • • 11 is the euclidean norm. This is then multiplied with a step function 
as given by equation (13), such that the force vanishes beyond some critical time. Alternatively one might 
find a solution going straight through a geometry with an angle of 7 t/ 3.5 using an upward force, 


F t = | y (°- 25 - ^ 2 ) 


Hr || < 0.5/ sin(7r/6) 
||r|j > 0.5/ sin(7r/6) 


5 Note that one has to set An ; i n i et = 1 when % 2 < e. 
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